clear
clear matrix
use "/.../ISSP 2006 RoG_Clean.dta", replace
set more off
xtset V3a


*recode and generate various variables needed for the analysis

gen unemspend=.
replace unemspend=1 if V23==1 | V23==2
replace unemspend=0 if V23==3 | V23==4 | V23==5

gen pensionspend=.
replace pensionspend=1 if V22==1 | V22==2
replace pensionspend=0 if V22==3 | V22==4 | V22==5


bysort V3a: egen mean_income=mean(income)
bysort V3a: gen rel_income=income/mean_income

gen rel_income_progpension=rel_income*progpension05
gen rel_income_stratunem=rel_income*stratunem05

gen pubsharepension=pubpension05/(pubpension05+privpension05)
gen rel_incpubsharepension=rel_income*pubsharepension
gen progpension_pubsharepension=progpension05*pubsharepension
gen rel_income_pubshareandprog=rel_income*pubsharepension*progpension05

gen pubsharepension_cat2=.
replace pubsharepension_cat2=1 if pubsharepension<.8
replace pubsharepension_cat2=2 if pubsharepension>.8 & pubsharepension<=1

gen empstat=.
*full-time
replace empstat=1 if WRKST==1
*part-time
replace empstat=2 if WRKST==2
*less than part-time, housewife/man, help family member, disabled, not in labor force
replace empstat=3 if WRKST==3 | WRKST==4 | WRKST==8 | WRKST==9 | WRKST==10
*unemployed
replace empstat=4 if WRKST==5
*in education
replace empstat=5 if WRKST==6
*retired
replace empstat=6 if WRKST==7






*regression models for Table 2 (comparison of pension and unemployment spending preferences)
xtset V3a

xi: xtlogit pensionspend rel_income female educ age i.empstat pubsharepension progpension05
outreg2 using "/.../Table2.doc", word replace

xi: xtlogit pensionspend rel_income female educ age i.empstat pubsharepension progpension05 rel_income_progpension
outreg2 using "/.../Table2.doc", word append

xi: xtlogit pensionspend rel_income female educ age i.empstat rel_income_progpension i.V3a 
outreg2 using "/.../Table2.doc", word append

xi: xtlogit unemspend rel_income female educ age i.empstat wfs stratunem05
outreg2 using "/.../Table2.doc", word append

xi: xtlogit unemspend rel_income female educ age i.empstat wfs stratunem05 rel_income_stratunem
outreg2 using "/.../Table2", word append

xi: xtlogit unemspend rel_income female educ age i.empstat rel_income_stratunem i.V3a 
outreg2 using "/.../Table2.doc", word append


*robustness check: models including past levels of spending (not shown in paper)

xi: xtlogit pensionspend rel_income female educ age i.empstat pubsharepension progpension05 rel_income_progpension
outreg2 using "/.../Table2a.doc", word replace

xi: xtlogit pensionspend rel_income female educ age i.empstat pubsharepension progpension05 pensions_macro rel_income_progpension
outreg2 using "/.../Table2a.doc", word append

xi: xtlogit unemspend rel_income female educ age i.empstat wfs stratunem05 rel_income_stratunem
outreg2 using "/.../Table2a.doc", word append

xi: xtlogit unemspend rel_income female educ age i.empstat wfs stratunem05 unempl_macro rel_income_stratunem
outreg2 using "/.../Table2a.doc", word append



*triple interaction, regression table (Table A.1 in the appendix)
xi: xtlogit pensionspend rel_income female educ age i.empstat progpension05 pubsharepension rel_income_progpension rel_incpubsharepension progpension_pubsharepension rel_income_pubshareandprog
outreg2 using "/.../TableA1.doc", word replace

xi: xtlogit pensionspend rel_income female educ age i.empstat rel_income_progpension i.V3a if pubsharepension_cat2==1
outreg2 using "/.../TableA1", word append

xi: xtlogit pensionspend rel_income female educ age i.empstat rel_income_progpension i.V3a if pubsharepension_cat2==2
outreg2 using "/.../TableA1", word append



*code for Figure 1, upper panel: interaction targeting and income for pensions 
*using the code supplied by Matt Golder: http://mattgolder.com/interactions#code

version 11.0
#delimit ;
set more off;

xi: xtlogit pensionspend rel_income progpension05 rel_income_progpension female educ age i.empstat pubsharepension;


matrix b=e(b);
matrix V=e(V);

scalar b1=b[1,1];
scalar b3=b[1,3];

scalar varb1=V[1,1];
scalar varb3=V[3,3];

scalar covb1b3=V[1,3];

scalar list b1 b3 varb1 varb3 covb1b3;

*     ****************************************************************  *;
*       Calculate data necessary for top marginal effect plot.          *;
*     ****************************************************************  *;

generate MVZ=_n;

replace  MVZ=. if _n>101;

gen conbx=b1+b3*MVZ if _n<101;

gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ) if _n<101;

gen ax=1.96*consx;

gen upperx=conbx+ax;

gen lowerx=conbx-ax;

*     ****************************************************************  *;
*       Construct stuff to produce the rug plot.  Need to create an     *;
*       offset position for the pipe marker, which will depend on       *;
*       where the histogram sits on the y-axis. This will requie some   *;
*       trial and error.                                                *;
*     ****************************************************************  *;

gen where=0.0;

gen pipe = "|";

egen tag_progpension = tag(progpension05);

*     ****************************************************************  *;
*       Construct variable to produce y=0 line.                         *;
*     ****************************************************************  *;

gen yline=0;

*     ****************************************************************  *;
*     ****************************************************************  *;
*       Produce marginal effect plot for X.                             *;
*     ****************************************************************  *;
*     ****************************************************************  *;

graph twoway hist progpension05, percent color(gs14) yaxis(2)
		||   scatter where progpension05 if tag_progpension, plotr(m(b 4)) ms(none) mlabcolor(gs5) mlabel(pipe) mlabpos(6) legend(off)
	    ||   line conbx   MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
        ||   line upperx  MVZ, clpattern(dash) clwidth(thin) clcolor(black)
        ||   line lowerx  MVZ, clpattern(dash) clwidth(thin) clcolor(black)
        ||   line yline  MVZ,  clwidth(thin) clcolor(black) clpattern(solid)
	    ||   ,
             xlabel(0 20 40 60 80 100, nogrid labsize(2))
		     ylabel(-0.8 -0.6 -0.4 -0.2 0.0, axis(1) nogrid labsize(2))
		     ylabel(0 10 20, axis(2) nogrid labsize(2))
	         yscale(noline alt)
		     yscale(noline alt axis(2))	
             xscale(noline)
             legend(off)
             xtitle("Pension progressivity" , size(2.5)  )
             ytitle("Marginal effect of income" , axis(1) size(2.5))
             xsca(titlegap(2))
             ysca(titlegap(2))
			 scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));

drop tag_progpension yline pipe where lowerx upperx ax consx conbx MVZ




*code for Figure 1, lower panel, interaction between targeting and income for unemployment 


version 11.0
#delimit ;
set more off;
*     ****************************************************************  *;
*     ****************************************************************  *;
*	  Get stuff ready to produce a marginal effect for X.               *;
*     ****************************************************************  *;
*     ****************************************************************  *;


xi: xtlogit unemspend rel_income stratunem05 rel_income_stratunem female educ age i.empstat wfs;

matrix b=e(b);
matrix V=e(V);

scalar b1=b[1,1];
scalar b3=b[1,3];

scalar varb1=V[1,1];
scalar varb3=V[3,3];

scalar covb1b3=V[1,3];

scalar list b1 b3 varb1 varb3 covb1b3;

*     ****************************************************************  *;
*       Calculate data necessary for top marginal effect plot.          *;
*     ****************************************************************  *;

generate MVZ=_n;

replace  MVZ=. if _n>51;

gen conbx=b1+b3*MVZ if _n<51;

gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ) if _n<51;

gen ax=1.96*consx;

gen upperx=conbx+ax;

gen lowerx=conbx-ax;

*     ****************************************************************  *;
*       Construct stuff to produce the rug plot.  Need to create an     *;
*       offset position for the pipe marker, which will depend on       *;
*       where the histogram sits on the y-axis. This will requie some   *;
*       trial and error.                                                *;
*     ****************************************************************  *;

gen where=0.0;

gen pipe = "|";

egen tag_stratunem = tag(stratunem05);

*     ****************************************************************  *;
*       Construct variable to produce y=0 line.                         *;
*     ****************************************************************  *;

gen yline=0;

*     ****************************************************************  *;
*     ****************************************************************  *;
*       Produce marginal effect plot for X.                             *;
*     ****************************************************************  *;
*     ****************************************************************  *;

graph twoway hist stratunem05, percent color(gs14) yaxis(2)
		||   scatter where stratunem05 if tag_stratunem, plotr(m(b 4)) ms(none) mlabcolor(gs5) mlabel(pipe) mlabpos(6) legend(off)
	    ||   line conbx   MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
        ||   line upperx  MVZ, clpattern(dash) clwidth(thin) clcolor(black)
        ||   line lowerx  MVZ, clpattern(dash) clwidth(thin) clcolor(black)
        ||   line yline  MVZ,  clwidth(thin) clcolor(black) clpattern(solid)
	    ||   ,
             xlabel(0 10 20 30 40 50, nogrid labsize(2))
		     ylabel(-0.8 -0.6 -0.4 -0.2 0.0, axis(1) nogrid labsize(2))
		     ylabel(0 10 20, axis(2) nogrid labsize(2))
	         yscale(noline alt)
		     yscale(noline alt axis(2))	
             xscale(noline)
             legend(off)
             xtitle("Low-income targeting of unemployment insurance" , size(2.5)  )
             ytitle("Marginal effect of income" , axis(1) size(2.5))
             xsca(titlegap(2))
             ysca(titlegap(2))
			 scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));

drop tag_stratunem yline pipe where lowerx upperx ax consx conbx MVZ







*Code for Figure 2: triple interaction between pension progressivity, public share and income


version 11.0
#delimit ;
set more off;
*     ****************************************************************  *;
*     ****************************************************************  *;
*	  Get stuff ready to produce a marginal effect for X.               *;
*     ****************************************************************  *;
*     ****************************************************************  *;


		
		
xi: xtlogit pensionspend rel_income progpension05 pubsharepension rel_income_progpension rel_incpubsharepension 
		progpension_pubsharepension rel_income_pubshareandprog
		female educ age;

		
		
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar b4=b[1,4];
scalar b5=b[1,5];
scalar b6=b[1,6];
scalar b7=b[1,7];

scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar varb4=V[4,4];
scalar varb5=V[5,5];
scalar varb6=V[6,6];
scalar varb7=V[7,7];

scalar covb1b4=V[1,4];
scalar covb1b5=V[1,5];
scalar covb1b7=V[1,7];
scalar covb4b5=V[4,5];
scalar covb4b7=V[4,7];
scalar covb5b7=V[5,7];

scalar list b1 b2 b3 b4 b5 b6 b7 varb1 varb2 varb3 varb4 varb5 varb6 varb7 
            covb1b4 covb1b5 covb1b7 covb4b5 covb4b7 covb5b7;

*     ****************************************************************  *;
*       Calculate data necessary for top marginal effect plot.          *;
*     ****************************************************************  *;

generate MVZ=(_n+45)/100;
replace MVZ=. if _n>55;

*     ****************************************************************  *;
*       We want to calculate the marginal effect of X on Y for all      *;
*       values of the modifying variable Z. We also want to calculate   *;
*       this marginal effect as Z changes for specific values of the    *;
*       second modifying variable W.  In the code below, we calculate   *;
*       the marginal effect of X on Y for all values of Z when W=0,     *;
*       when W=1, when W=2, and when W=3.                               *;
*     ****************************************************************  *;

#delimit ;
gen conb0=b1+b4*0+b5*MVZ+b7*(0*MVZ) if _n<55;
gen conb1=b1+b4*25+b5*MVZ+b7*(25*MVZ) if _n<55;
gen conb2=b1+b4*50+b5*MVZ+b7*(50*MVZ) if _n<55;
gen conb3=b1+b4*75+b5*MVZ+b7*(75*MVZ) if _n<55;
gen conb4=b1+b4*100+b5*MVZ+b7*(100*MVZ) if _n<55;

*     ****************************************************************  *;
*       Calculate the standard errors for the marginal effect of X on Y *;
*     ****************************************************************  *;

gen conse0=sqrt(varb1
                + varb4*(0^2) + varb5*MVZ^2 + varb7*MVZ^2*(0^2)
                + 2*0*covb1b4 + 2*MVZ*covb1b5 + 2*0*MVZ*covb1b7+2*0*MVZ*covb4b5
                + 2*(0^2)*MVZ*covb4b7 + 2*0*(MVZ^2)*covb5b7)  if _n<55;
                
gen conse1=sqrt(varb1
                + varb4*(25^2) + varb5*MVZ^2 + varb7*MVZ^2*(25^2)
                + 2*25*covb1b4 + 2*MVZ*covb1b5 + 2*25*MVZ*covb1b7+2*25*MVZ*covb4b5
                + 2*(25^2)*MVZ*covb4b7 + 2*25*(MVZ^2)*covb5b7)  if _n<55;
                
gen conse2=sqrt(varb1
                + varb4*(50^2) + varb5*MVZ^2 + varb7*MVZ^2*(50^2)
                + 2*50*covb1b4 + 2*MVZ*covb1b5 + 2*50*MVZ*covb1b7+2*50*MVZ*covb4b5
                + 2*(50^2)*MVZ*covb4b7 + 2*50*(MVZ^2)*covb5b7)  if _n<55;
                
gen conse3=sqrt(varb1
                + varb4*(75^2) + varb5*MVZ^2 + varb7*MVZ^2*(75^2)
                + 2*75*covb1b4 + 2*MVZ*covb1b5 + 2*75*MVZ*covb1b7+2*75*MVZ*covb4b5
                + 2*(75^2)*MVZ*covb4b7 + 2*75*(MVZ^2)*covb5b7)  if _n<55;
                
gen conse4=sqrt(varb1
                + varb4*(100^2) + varb5*MVZ^2 + varb7*MVZ^2*(100^2)
                + 2*100*covb1b4 + 2*MVZ*covb1b5 + 2*100*MVZ*covb1b7+2*100*MVZ*covb4b5
                + 2*(100^2)*MVZ*covb4b7 + 2*100*(MVZ^2)*covb5b7)  if _n<55;
                
*     ****************************************************************  *;
*                           Create t statistics                         *;
*     ****************************************************************  *;

gen t0=conb0/conse0;
gen t1=conb1/conse1;
gen t2=conb2/conse2;
gen t3=conb3/conse3;
gen t4=conb4/conse4;

*     ****************************************************************  *;
*       Generate a `shadow' variable that is equal to the marginal      *;
*       effects that have already been calculated.  These will be       *;
*       helpful in placing the stars in the figure as you will see.     *;
*     ****************************************************************  *;

gen consb0=conb0;
gen consb1=conb1;
gen consb2=conb2;
gen consb3=conb3;
gen consb4=conb4;

*     ****************************************************************  *;
*       Make the `shadow' variable be missing if the t score is not     *;
*       larger than the critical level of significance that you want.   *;
*     ****************************************************************  *;

replace consb0 = . if abs(t0)<2.01;
replace consb1 = . if abs(t1)<2.01;
replace consb2 = . if abs(t2)<2.01;
replace consb3 = . if abs(t3)<2.01;
replace consb4 = . if abs(t4)<2.01;

*     ****************************************************************  *;
*       Generate a string variable called str1 that is designated with  *;
*       a star.                                                         *;
*     ****************************************************************  *;

generate str1 txt="*";

set textsize 100;

graph twoway line conb0 MVZ, clpattern(solid) clwidth(thin) 
        || scatter consb0 MVZ, mlabel(txt) msymbol(i) mlabsize(vsmall) mlabgap(-1.0) mlabposition(11)  
        || line conb1 MVZ, clpattern(solid) clwidth(thin) lcolor (blue)
        || scatter consb1 MVZ, mlabel(txt) msymbol(i) mlabsize(vsmall) mlabgap(-1.0) mlabposition(11)  
        || line conb2 MVZ, clpattern(solid) clwidth(thin) lcolor (red)
        || scatter consb2 MVZ, mlabel(txt) msymbol(i) mlabsize(vsmall) mlabgap(-1.0) mlabposition(11) 
        || line conb3 MVZ, clpattern(solid) clwidth(thin)  lcolor (mint)
        || scatter consb3 MVZ, mlabel(txt) msymbol(i) mlabsize(vsmall) mlabgap(-1.0) mlabposition(11)          
        || ,       
        xlabel(.45 .50 .55 .60 .65 .70 .75 .80 .85 .90 .95 1, nogrid labsize(2.5))   
        ylabel(-.8 -.6 -.4 -.2 0 .2, nogrid labsize(2.5))  
        yscale(noline)  xscale(noline)  legend(off)  
        yline(0, lcolor(black))   
        xtitle("Public Share in Pension Spending")  
        ytitle("Marginal Effect of Income")  
        xsca(titlegap(4)) ysca(titlegap(4)) 
        scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));

		
drop MVZ conb0 conb1 conb2 conb3 conb4 conse0 conse1 conse2 conse3 conse4 t0 t1 t2 t3 t4 consb0 consb1 consb2 consb3 consb4 txt

log close;



